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Abstract 

Traditional  transmission  travel-time  tomography  hinges  on  ray  tracing  techniques. 
We  propose  a  PDE-based  Eulerian  approach  to  travel-time  tomography  so  that  we 
can  avoid  the  cumbersome  ray-tracing.  We  start  from  the  eikonal  equation,  define 
a  mismatching  functional  and  derive  the  gradient  of  the  nonlinear  functional  by  an 
adjoint  state  method.  The  resulting  forward  and  adjoint  problems  can  be  efficiently 
solved  by  using  the  fast  sweeping  method;  a  limited  memory  BFGS  method  is  used 
to  drive  the  mismatching  functional  to  zero  with  quadratic  convergence.  2-D,  3-D 
numerical  results  and  Marmousi  synthetic  velocity  model  demonstrate  the  robustness 
and  the  accuracy  of  the  method. 


1  Introduction 

Estimation  of  wave-speed  distribution  from  acoustic,  seismic  or  electromagnetic  first-arrival 
travel-time  data  is  the  goal  of  transmission  travel-time  tomography.  In  seismics  velocity 
analysis  is  often  an  important  step  in  prospect  evaluation  in  areas  where  lithology  and 
structure  undergo  significant  lateral  change.  In  this  work  we  propose  a  new,  robust  and 
efficient  tomography  method  which  is  aimed  at  such  applications. 

All  the  traditional  methods  of  travel-time  tomography  are  directly  based  on  Fermat’s 
least  travel-time  principle  and  bear  a  close  link  to  the  X-ray  computerized  tomography  (CT) 
used  in  medical  diagnosis.  In  medical  CT  the  measured  data  are  assumed  to  be  modeled 
by  line  integrals  of  wave  amplitude  attenuation  for  straight  ray-paths  passing  through  the 
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body,  and  the  Radon  transform  provides  the  foundation  for  medical  CT.  However,  in  seismics 
the  ray-path  curvature  has  to  be  taken  into  account  in  that  lithology  and  structure  usually 
have  strong  inhomogeneity,  and  the  resulting  ray-paths  can  depend  strongly  on  the  unknown 
wave  speeds.  To  achieve  such  a  purpose,  ray-tracing  based  travel-time  tomography  methods 
require  very  complicated  data  structure  to  trace  curved  rays  through  each  pixel  [4];  see  [25] 
for  3-D  examples.  In  addition,  such  ray-tracing  based  methods  inevitably  produce  irregular 
ray  coverage  of  the  computational  domain,  and  the  resulting  system  of  equations  may  not 
be  well-conditioned  [1,  2,  3].  In  this  paper  we  propose  a  PDE-based  Eulerian  approach  to 
travel-time  tomography  so  that  we  can  avoid  the  cumbersome  ray-tracing. 

Recall  that  a  necessary  condition  for  Fermat’s  least  travel-time  principle  to  hold  is  charac¬ 
terized  by  the  eikonal  equation  for  travel-time  [11],  and  the  viscosity  solution  for  the  eikonal 
equation  with  a  point-source  condition  is  the  least  travel-time  from  the  source  to  an  arbitrary 
point  connected  by  a  shortest  ray-path,  as  observed  by  [13,  16].  Because  of  its  continuous  de¬ 
pendence  on  the  wave-speed  distribution  and  source  locations,  the  viscosity  solution  can  be 
computed  by  various  numerical  schemes  stably.  In  this  work,  we  model  travel-times  from  a 
single  source  to  multiple  receivers  by  using  the  eikonal  equation  and  propose  a  fast  sweeping 
based  adjoint  state  method  for  transmission  travel-time  tomography.  The  new  approach  not 
only  overcomes  some  shortcomings  inherited  in  the  traditional  ray-tracing  based  travel-time 
tomography  but  also  enjoys  quadratic  convergence,  thus  it  is  very  fast  and  robust. 

However,  first-arrival  based  transmission  travel-time  tomography  usually  has  very  lim¬ 
ited  resolution.  Since  the  output  from  travel-time  tomography  is  mainly  used  for  building 
a  macro  velocity  model  in  seismic  velocity  analysis,  it  is  important  to  have  very  fast,  effi¬ 
cient  tomography  tools,  even  though  we  may  have  to  use  only  first-arrivals.  On  the  other 
hand,  multi-valued  travel-times  and  resulting  multipathings  are  common  in  complex  veloc¬ 
ity  structures,  it  is  necessary  to  take  into  account  all  the  arrivals  systematically.  As  is  well 
known,  ray-tracing  methods  can  yield  all  arrivals,  and  the  works  presented  in  [8,  7]  have  used 
multivalued  travel-times  from  ray  tracing  methods,  but  those  works  are  on  reflection  travel¬ 
time  tomography  which  is  different  from  transmission  tomography  in  that  rays  start  at  the 
surface,  reflect  off  interfaces  whose  depths  are  to  be  determined,  and  return  to  the  surface. 
To  use  all  arrivals  in  transmission  tomography  in  an  Eulerian  framework,  we  proposed  to 
formulate  transmission  tomography  by  using  the  Liouville  equation  based  PDE  framework 
in  phase  space  [15,  14], 

In  this  work,  we  will  only  concentrate  on  first-arrival  based  traveltime  tomography.  We 
start  from  a  mismatching  functional  between  measured  and  simulated  data  and  drive  the 
functional  to  zero  by  a  well  designed  limited  memory  BFGS  (L-BFGS)  optimization  method. 
Although  our  approach  shares  some  similarities  with  that  in  Sei-Symes  [21,  22],  that  work  was 
based  on  the  paraxial  formulation  of  the  eikonal  equation  and  only  illustrated  the  feasibility  of 
computing  the  travel-time  gradient  by  using  the  adjoint  state  method.  Instead  of  the  paraxial 
formulation  of  the  eikonal  equation,  we  derive  the  gradient  of  the  mismatching  functional 
directly  from  the  steady  eikonal  equation  by  using  the  adjoint  state  method;  furthermore, 
we  apply  the  fast  sweeping  method  [26,  24,  12]  to  solve  the  eikonal  equation  directly  (the 
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forward  problem)  and  design  a  new  fast  sweeping  method  to  solve  the  adjoint  equation  of 
the  linearized  eikonal  equation  (the  adjoint  problem)  so  that  the  required  gradient  can  be 
computed  highly  efficiently;  finally  a  limited  memory  BFGS  optimization  method  drives  the 
mismatching  functional  to  zero  with  quadratic  convergence. 

In  Section  2,  we  first  explain  the  variational  method  for  inverting  the  velocity  model 
using  measurements  on  the  boundary  of  a  specified  computational  domain.  To  minimize  the 
energy  in  the  variational  formulation,  we  derive  the  gradient  of  the  nonlinear  functional.  To 
efficiently  compute  the  gradient  direction,  in  Section  3  we  apply  the  fast  sweeping  method 
to  the  eikonal  equation  and  design  a  new  fast  sweeping  method  for  the  adjoint  equation  of 
the  linearized  eikonal  equation.  Sections  4,  5  and  6  show  various  numerical  examples  to 
demonstrate  the  feasibility  and  the  robustness  of  the  new  formulation.  Section  7  will  then 
conclude  the  paper. 

2  Governing  Equations 

We  start  from  the  eikonal  equation  with  a  point  source  condition  in  an  isotropic  medium 
which  occupies  an  open,  bounded  rectangular  domain  Qp  C  TZ3.  By  isotropy  here  we  mean 
the  wave  velocity  has  no  directional  dependence.  The  equation  is  as  follows, 

|VT|  =  I  (1) 

with  the  point  source  condition 

T(xs)  =  0,  (2) 

where  T(x)  is  the  travel-time  of  wave  from  the  source  xs  to  the  point  x,  and  c  G  C'1(f2)  is  a 
positive  velocity  function. 

For  a  given  velocity  model  c,  the  viscosity  solution  of  this  equation  can  be  computed 
efficiently  by  fast  sweeping  methods,  and  such  solutions  correspond  to  the  least  travel-time 
or  the  first-arrival  travel-time  according  to  [16]. 

In  this  work  we  are  interested  in  the  related  inverse  problem,  the  so-called  transmission 
travel-time  tomography  problem:  given  both  the  first-arrival  travel-time  measurements  on 
the  boundary  dfl.p  and  the  location  of  the  point  source  xs  G  Qp,  invert  for  the  velocity  held 
c(x)  inside  the  domain  Qp. 

To  achieve  this  we  propose  to  invert  for  the  velocity  model  by  minimizing  the  following 
mismatching  functional  (energy), 

E(c)  =  \j  \T-rf-,  (3) 

"  J  dftp 

where  T*\qqp  is  the  measurement  and  T\qqp  is  computed  by  solving  (1)  with  a  point  source 
condition  (2).  In  other  words,  this  energy  measures  the  L2- difference  between  the  experi¬ 
mental  measurement,  T*,  and  the  solution  from  the  eikonal  equation,  T,  on  the  boundary 
of  the  computational  domain. 
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To  minimize  this  energy,  we  use  the  method  of  gradient  descent.  We  first  perturb  the 
velocity  field  c  by  ec,  which  causes  a  corresponding  change  in  T  by  eT.  The  change  in  the 
energy  is  then  given  by 

5E  =  e  [  T(T  —  T*)  +  0(e2) .  (4) 

Jactp 

From  the  state  equation  (1),  the  perturbations  in  c  and  T  are  related  by 

TXTX  +  TyTy  +  TZTZ  =  —  —  .  (5) 

We  need  to  determine  the  perturbation  in  c,  c,  so  as  to  decrease  the  energy  E(c).  The 
main  difficulty  is  that  the  perturbation  in  E,  5E,  depends  on  c  implicitly  through  T  and  the 
partial  differential  equation  (5).  To  efficiently  compute  c  which  minimizes  E,  we  apply  the 
adjoint  state  method. 

Multiplying  (5)  by  eA,  integrating  it  over  ttp,  applying  integration  by  parts,  and  adding 
the  resulting  expression  to  (4),  we  have 

—  =  /  T(T-T>)+  f  j  +  f  f  \T,T\l=  +  f  f  \rt\l— 

^  J  d£lp  J  y  J  z  J  x  J  z  J  x  J  y 

-  f  T[(XTx)x  +  (A Ty)y  +  (A Tz)z]  +  f  —  +  0(e) .  (6) 

Jflp  Jflp  c 

Next,  we  choose  A  satisfying 


{~Tx)X\x  +  [(— Ty)X]y  +  [(— TZ)X]Z  —  0  , 


(7) 


with  the  boundary  condition, 

(n  •  VT)  A  =  T*-T,  (8) 

on  the  boundary  d Qp,  where  n  is  the  unit  outward  normal  of  the  boundary.  By  introducing 
this  adjoint  state  equation,  one  can  eliminate  the  dependence  of  T  when  determining  the 
gradient  of  E  with  respect  to  c. 

Ignoring  the  higher  than  linear  order  terms  in  the  energy  perturbation,  we  have 

:  //'• 


To  minimize  the  energy  using  the  method  of  gradient  descent,  one  could  choose  the 
perturbation  c—  —  A/c3.  This  implies 


SE  —  —e  [  c2  <  0  (10) 

J  Qp 

and  the  equality  holds  when  ||c||#o(n  )  =  0.  However,  it  is  not  straight-forward  how  one  can 
guarantee  the  following  two  properties, 
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1-  ck\dnv  =  0; 

2.  ck+1  =  ck  +  eck  smooth. 

The  first  condition  assumes  that  we  can  measure  c  on  the  boundary  dflp,  denoted  by 
c*\dnp >  which  is  a  reasonable  assumption.  This  means  that  the  variations  of  the  velocity 
function  on  the  boundary  should  be  zero. 

The  second  condition  is  a  regularity  condition  on  ck.  This  regularity  seems  to  be  too 
restrictive  in  practice.  In  general,  one  only  needs  ck  G  C1  to  guarantee  well-posedness  of 
the  state  equation  (1).  However,  assuming  that  one  uses  ck  =  — A/c3  directly,  it  is  not  clear 
whether  this  function  would  give  us  the  desired  regularity.  Even  if  this  perturbation  is  in  C1i 
the  numerical  solution  may  have  jumps  or  spikes.  These  irregularities  will  force  one  to  pick  a 
very  small  step-size,  ek,  in  the  minimization  process.  Therefore,  to  have  faster  convergence, 
we  impose  the  above  regularity  in  each  iteration. 

One  way  to  satisfy  the  above  two  properties  is  to  use  the  descent  direction 

c  =  -(/  -  ^A)-1  (A)  ,  dl) 

where  I  is  the  identity  operator,  A  is  the  Laplacian  operator  and  v  >  0  controls  the  amount 
of  regularity  that  one  wants.  The  homogeneous  boundary  condition  is  imposed  in  inverting 
the  operator  (/  —  z/A).  With  this  particular  c,  we  have 


8E  —  —e  f  (c2  +  z/|Vc|2)  <  0.  (12) 

Jnp 

We  notice  that  this  process  amounts  to  seeking  updates  in  some  weighted  Sobolev  space  in 
the  case  v  >  0.  Then  the  above  equality  holds  when  ||c||//i(n  )  =  0. 

In  the  above  calculation,  we  use  the  first- arrivals  at  different  receivers  associated  with  a 
single  point  source.  If  we  perform  multiple  such  experiments,  namely,  we  have  many  such 
data  sets,  then  those  can  be  easily  incorporated  into  the  formulation.  For  example,  we  can 
assume  that  there  are  N  point  sources  located  at  x®,  for  i  =  1,  •  •  • ,  N,  and  N  sets  of  first- 
arrival  travel-time  measurements  T*  associated  with  these  N  sources  are  available.  Then  we 
can  simply  define  a  new  energy 

1  N  f 

=  A  /  <13) 

Z  i= l 


where  Tt  is  the  solution  from  the  eikonal  equation  with  the  corresponding  point  source  con¬ 
dition  T(x()  =  0.  Utilizing  the  same  approach  as  above,  we  have  the  following  perturbation 
in  the  energy 


8EN 


_  N 


i=  1 


6 


(14) 
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where  A  $  is  the  adjoint  variable  of  T*  satisfying 

{  [ -  (2^)3;]  Ai}a;  +  {[—  {Ti)y\\i]y  +  {  [~  ( Tj  )  3]  Aj  }  3  =  0  ,  (l5) 

with  the  boundary  condition, 

(n  •  VTj)Aj  =  T*  —  Tj ,  (16) 

for  i  —  1,  •  •  • ,  N. 

Consequently,  we  can  choose  the  following  gradient  direction  to  minimize  the  energy 
EN(c), 

c  =  -(/-„A)-dif>V  (17) 

We  remark  that  the  above  updating  procedure  is  similar  to  the  so-called  simultaneous 
iterative  reconstruction  technique  frequently  used  in  medical  imaging;  it  is  also  possible  to 
adopt  the  algebraic  reconstruction  type  technique  as  used  in  [9]  to  update  the  velocity. 

3  Algorithm  and  Numerical  Implementations 

3.1  Tomography  algorithm 

Here  we  give  an  algorithm  for  this  tomography  problem. 

Tomography  Algorithm: 

1.  Initialize  ck  for  k  =  0  by  solving 

(/  -  v A )c°  =  0  ,  (18) 

with  the  boundary  condition  c°|,9q  =  cexact|,gn- 

2.  Obtain  T(x,z )  by  solving  (1)  with  the  point  source  condition  (2)  using  c—  ck. 

3.  Obtain  A (x,z)  by  solving  (7)  with  the  boundary  condition  (8). 

4.  Obtain  ck  using  (11). 

5.  Determine  ek  using,  for  example,  the  Armijo-Goldstein  rule  or  simply  ek  =  e. 

6.  Update 

ck+i  =  ck  +  ek£k  _ 

7.  Go  back  to  Step  2  until  \\ck{x,  z)\\2  <  5  or  k  >  kmax,  where  <5  and  kmax  are  given 
convergence  parameters. 
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To  start  the  iteration,  we  need  to  initialize  c°.  Here  in  the  algorithm,  we  assume  that 
we  can  measure  the  velocity  at  receivers,  giving  c°|an  =  cexact  I  on  ■  This  condition  can 
be  replaced  by  other  assumptions.  In  practice,  due  to  the  nonlinearity  of  the  problem, 
different  initial  guesses  will  generally  converge  to  different  energy  minimizers.  This  non¬ 
uniqueness  can  be  overcome  by  some  a  priori  knowledge  of  the  model.  For  example,  the 
above  assumption  can  be  relaxed  by  replacing  the  Dirichlet  conditions  on  both  the  left  and 
right  boundaries  with  the  Neumann  boundary  conditions  dc°  /  dx\x=Xmin  =  dc° / dx\x=Xm^  =  0. 

The  convergence  could  be  sped  up  by  replacing  the  gradient  descent  method  with  BFGS- 
type  iterations.  To  solve  the  elliptic  equation  in  Step  4,  we  use  the  FFT.  In  Step  2  and  Step 
3,  both  the  equations  (1)  and  (7)  can  be  solved  by  the  Fast  Sweeping  Method  [26,  24,  12], 
which  we  detail  next. 


3.2  Fast  Sweeping  Method  for  Equation  (1) 

The  fast  sweeping  method  was  originated  in  Boue  and  Dupis  [5],  its  first  PDE  formulation 
was  in  implicit  and  non-parametric  shape  reconstruction  from  unorganized  points  using  a 
variational  level  set  method  [27];  Zhao  [26]  proved  the  O(N)  convergence  of  the  method 
for  the  eikonal  equation  based  on  the  Godunov  Hamiltonian  on  Cartesian  meshes;  later 
on,  the  fast  sweeping  method  was  extended  to  treat  Hamilton- Jacobi  equations  with  convex 
Hamiltonians  based  on  the  Godunov  Hamiltonian  [24]  and  handle  Hamilton-Jacobi  equations 
with  non-convex  Hamiltonians  based  on  the  Lax-Friedrichs  Hamiltonian  [12];  see  [24,  12]  and 
references  therein  for  the  fast  sweeping  method  on  Cartesian  meshes  and  [20]  for  the  method 
on  triangulated  meshes.  Certainly,  one  may  also  use  other  methods  such  as  the  fast  marching 
method  [23]. 

To  be  self-contained,  we  give  a  short  summary  of  the  fast  sweeping  method  for  eikonal 
equations.  To  avoid  cluttered  notations  we  present  the  algorithm  for  the  2-D  case  only;  see 
[26]  for  more  details. 

First  we  discretize  the  rectangular  domain  C  77  2  into  a  uniform  mesh  with  mesh  points 
xy  and  mesh  sizes  Ax  =  Az  =  h,  and  we  denote  the  numerical  solution  at  Xjj  by  TJj. 
Applying  the  Godunov  numerical  Hamiltonian  to  the  eikonal  equation,  for  i  —  2,  —  1, 

j  —  2,  •  •  • ,  J  —  1,  we  have 


[(Tij  -  Txmin)+r  +  -  Tzmin)+Y  = 


h 2 


h3 


(19) 


where 


'-h'xmin  min(Tj_lj-,  Tj+lj),  Tzmin  min(T]  j_j,  TJj-i-i), 

and  (a;)+  denotes  the  positive  part  of  x.  At  the  boundary  of  the  computational  domain  one 
sided  difference  is  used. 


Fast  Sweeping  Algorithm: 
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1.  Initialize  the  point  source  condition  T(xs)  =  0  by  assigning  the  exact  value  if  xs  is  a  mesh 
point,  or  assigning  to  grid  points  near  xs  exact  values  which  are  computed  by  using  the 
constant  velocity  at  the  point  source.  These  values  are  fixed  in  later  iterations.  Assign 
larger  positive  values  at  all  other  grid  points,  and  these  values  will  be  updated  later. 

2.  Update  the  solution  by  Gauss-Seidel  iterations  with  alternating  sweeping.  At  each  grid  point 

x,j  whose  value  was  not  fixed  during  the  initialization,  compute  the  candidate  solution, 
denoted  by  T  of  (19)  from  the  current  values  of  its  neighbors  Ttj±l  and  then  update 

Tij  to  be  the  smaller  one  between  T  and  its  current  value;  i.e. ,  T™w  =  niin(T°jd,  T).  We 
sweep  the  whole  domain  with  four  alternate  ordering  repeatedly:  i  —  1  :  I,j  —  1  :  J; 
i  —  1  :  I,  j  —  J  :  1;  i  —  I  :  1,  j  —  1  :  J;  i  —  I  :  1,  j  —  J  :  1.  Here  i  and  j  are  the  running 
indices  along  x  and  y  directions. 

3.  Test  the  convergence:  given  convergence  criterion  e  >  0,  check  whether  ||Tn+1  —  Tn\\Li  < 
e. 

We  remark  that  the  sweeping  strategy  can  be  used  for  more  general  Hamilton- Jacobi 
equations  as  long  as  an  efficient  local  solver  is  available  at  each  grid  point,  so  that  an  iterative 
procedure  is  well  defined  at  each  local  grid  point.  On  the  other  hand,  we  may  apply  the 
sweeping  strategy  to  solve  equation  (7),  which  reduces  to  a  symmetrical  Gauss-Seidel- type 
iteration  method  . 


3.3  Fast  Sweeping  Method  for  Equation  (7) 

Next  we  design  a  fast  sweeping  method  for  equation  (7).  Once  again  to  simplify  the  notation, 
we  give  a  2-D  formulation  only;  the  extension  to  a  3-D  formulation  is  straightforward. 

The  adjoint  state  equation  (7)  can  be  written  in  the  following  form 

(a\)x  +  (b\)z  =  0  ,  (20) 


where  a  and  b  are  given  functions  of  (x,  z ). 

Considering  a  computational  cell  centered  at  (xtl  Zj)  and  discretizing  the  equation  in 
conservation  form,  we  have 


(<h+l/2jAj+i/2j  —  Oj_i/2jAj_i/2j)  +  (bij+i/2\i,j+l/2  ~  Ajj_i/2)  —  0  .  (21) 


The  values  of  A  on  the  interfaces,  Xi±i/2,j  and  \ij±i/2,  are  determined  according  to  the  prop¬ 
agation  of  characteristics.  In  the  case  when  o*+i/2j  >  0,  the  characteristic  for  determining  A 
goes  from  the  left  hand  side  of  the  interface  to  the  right  hand  side,  and  this  suggests  that  we 
use  the  value  A*j  to  define  \t+i/2,y,  otherwise,  we  have  A,+i/2j  =  Aj+ij.  The  terms  Xij±i/2 
can  be  defined  in  a  similar  way. 
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Introducing  the  following  notations 

n±  _  «i+l/2J  =t  |«j+l/2j| 
ai+l/2J  ~  2 

,  ±  _  Kj+1/2  ±  \bjj+i/2\ 

°hj+ 1/2  _  2 

we  have 


^7  (iK,j+l/2Xtj  ~~  \j+l/2Xi,j+1)  —  (K,j+l/2Xtj-l  —  \j+l/2Xi,j)  )  =  0  >  (22) 

which  can  be  rewritten  as 


ai-l/2,j  i  |Oj-l/2j| 


A  _ _ 

h— 1/2,/  2 


and  btj~  1/2 


bi,j- 1/2  =t  1/2  | 


ai+l/2,j  ai-  1/2J  ^J+l/2  bij~\/2 


Ax 


A  z 


Xi,j  ~ 


ai-l/2,jXi~1d  ai+l/2,iX 


i+l/2,jAi+1J 


Ax 

|  ^tj-l/2^0^1  ~  ftj,j+l/2-^J+l  /23> 

Ax; 


This  gives  an  expression  to  build  up  a  fast  sweeping-type  iterative  method. 

To  apply  this  iterative  scheme  to  equation  (7),  we  need  to  specify  the  function  values  of 
a  and  b  not  at  the  cell  centers  (ay,  Zj),  but  on  the  cell  interfaces  (ay±1/2 ,Zj)  and  (xi,Zj: tl/2). 
This  can  be  done  easily  using  central  differences.  For  example,  we  have  a*+i/2j  =  —(Ti+ij  — 
Titj)/ Ax  and  at-i/2,j  =  —  Tj_ij)/Ax.  In  addition,  we  have  to  incorporate  the  boundary 

condition  (8)  into  the  above  linear  system  for  A  as  well.  Then  we  can  show  that  the  coefficient 
matrix  of  the  resulting  linear  system  for  A  is  irreducibly  diagonally  dominant,  therefore  the 
alternating  symmetrical  Gauss-Seidel  iteration  converges. 


Fast  Sweeping  Algorithm  for  equations  (7)  and  (8): 

1.  On  the  boundary,  compute  (n  ■  VT)  from  the  solution  of  the  eikonal  solver  using  one  side 
difference.  Next,  compute  the  boundary  condition  for  A  according  to  (8).  These  values  will 
be  fixed  in  the  following  computations. 

2.  Update  \l)3  at  the  interior  points  according  to  (22).  As  in  the  Fast  Sweeping  method  for 
(1),  we  sweep  the  whole  domain  with  four  alternate  orderings. 

3.  For  some  given  convergence  criterion  e  >  0,  repeat  2  until  ||An+1  —  A"||li  <  e. 

We  point  out  that  the  above  fast  sweeping  method  is  different  from  the  fast  marching 
method  used  in  [10],  in  that  our  method  is  iterative  and  theirs  is  constructive  based  on 
upwinding  properties. 
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3.4  L-BFGS  method 

In  the  Tomography  Algorithm,  we  update  the  approximation  to  the  velocity  by  the  typical 
gradient  descent  method,  where 

ck+i  =  ck  _  ek~k  _  (24) 

Although  it  is  simple  to  implement,  the  method  is  not  efficient  because  it  takes  a  large 
number  of  iterations  to  converge  to  the  steady  state  solution. 

To  speed  up  the  convergence,  we  can  apply  the  quasi-Newton  method  defined  by 

ck+ 1  =  ck  +  eksk  ,  (25) 

where  sk  =  —A^lE\ck)  and  Ak  is  a  positive  definite  operator  satisfying  the  secant  condition 

Ak+1(ck+1  -  ck)  =  E\ck+l )  -  E'(ck ) .  (26) 

In  this  iteration,  the  operator  Ak+\  is  updated  by  modifying  the  previous  operator  Ak. 

One  possible  way  to  modify  this  operator  is  defined  by  Broydon-Fletcher-Goldfarb- 
Shanno  (BFGS)  procedure, 

Akv  =  Ak_iv  +  a  <  p,v  >  p  +  (3  <  q,v  >  q ,  (27) 

where 

p  =  y/\\y\\  ,  q  =  Ak-is/\ \ Ak_xs\ \ , 

a  =  \\y\\2/ <  y,s  >  ,  (3  =  -| \Ak_is\  j J/  <  s,  Ak_is  >  (28) 

with  s  =  ck  —  cfc_1,  y  =  E'(ck)  —  A'(cfc_1)  and  A0  =  L 

However,  in  practice,  the  condition  number  of  Ak  can  be  increased  significantly  through¬ 
out  the  iteration,  which  makes  the  computation  inaccurate.  To  alleviate  this,  one  can  modify 
the  iteration  using  the  limited  memory  BFGS  (L-BFGS)  given  by 

k 

Akv  =  v+  ^2  (aj  < Pj> v  > Pj  +  Pj  <  qp v  >  qj)  ■  (29) 

j=k—L+ 1 

In  this  paper,  we  adapt  the  L-BFGS-B  code  from  [6].  This  code  requires  user  to  provide  only 
subroutines  to  compute  both  the  energy  to  be  minimized  and  the  gradient  of  this  energy. 
The  step  size  ek  is  automatically  determined. 

4  Two-Dimensional  Numerical  Examples 

In  the  following  examples,  we  use  129  x  129  grid  points  in  the  x-z  space.  Using  the  above 
formulation,  we  need  measurements,  denoted  by  <j>*  and  T*,  on  the  boundary  diip.  If  the 
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point  source  is  located  inside  fl,  the  characteristics  of  the  eikonal  equation  always  flow  out 
from  the  domain.  Therefore,  in  synthetic  experiment  the  boundary  measurements  can  be 
obtained  by  solving  the  equation  (1)  directly  using  the  Fast  Sweeping  Method  together  with 
the  exact  velocity  c. 

For  each  velocity  model  below,  we  have  implemented  the  following  two  cases  -  one  source 
and  ten  sources.  For  the  one  source  case,  we  use  the  boundary  measurements  from  the 
only  point  source  located  at  (x,z)  =  (0,0.1).  In  the  cases  with  ten  point  sources,  we  use 
nine  more  sets  of  boundary  measurements,  and  these  correspond  to  source  locations  at 
(x,  z)  =  (±0.25,  0.1),  (±0.5,  0.1),  (0, 1.9),  (±0.25, 1.9)  and  (±0.5, 1.9),  respectively.  However, 
to  save  some  space  we  only  present  the  results  corresponding  to  the  case  of  ten  sources. 

To  start  the  algorithm,  we  initialize  the  velocity  c°  by  solving  the  above  elliptic  equation 
(18)  with  v  —  1. 

4.1  Example  1.  Constant  Model. 

The  exact  velocity  model  is  given  by  c  =  1. 

We  use  the  BFGS  method  to  invert  for  the  velocity.  The  results  are  shown  in  Figure  1.  As 
we  can  see,  the  recovered  velocity  is  almost  exact,  the  relative  error  is  almost  negligible,  and 
we  observe  the  typical  quadratic  convergence  of  the  algorithm  due  to  the  L-BFGS  method. 

4.2  Example  2.  Waveguide  Model. 

The  exact  velocity  model  is  given  by 

c(x,  z)  =  3  —  2.5  exp  ^  .  (30) 

We  apply  the  BFGS  method  to  invert  for  the  velocity.  Figure  2  shows  the  relative  error 
at  convergence  and  the  convergence  history.  In  Figure  3,  we  show  slices  of  the  cross-sections 
of  the  solution  along  z  =  1  and  x  =  0.  As  we  can  see,  we  are  able  to  recover  the  velocity  as 
well. 


4.3  Example  3.  Gaussian  Model. 

The  exact  velocity  model  is  given  by 


c(x ,  z)  —  3  —  -  exp 


x2  ±  (z  —  0.5)" 
052 


—  exp 


x2  +  (z- 1.25  y 
052 


(31) 


Next  we  apply  the  BFGS  method  to  invert  for  the  velocity.  In  Figure  4,  we  show  the 
convergent  velocity  and  the  convergence  history  of  the  algorithm;  once  again,  we  observe 
quadratic  convergence.  In  Figure  5,  we  show  slices  of  cross-sections  of  the  final  converged 
velocity;  as  we  can  see,  they  fit  well  with  the  exact  velocity. 
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(a) 


(b) 


(c) 


(d) 


Figure  1:  (Example  1.  Ten  Sources)  BFGS.  (a):  the  initial  guess;  (b):  final  approximated  c; 
(c):  the  relative  error  in  the  solution;  (d):  the  convergence  history  of  energy. 


To  further  test  the  algorithm,  we  repeat  the  experiment  but  perturb  the  synthetic  data 
T*  with  some  noise.  Using  the  same  velocity  model,  we  first  compute  the  travel-time  on 
the  boundary  of  the  domain.  These  measurements  are  added  5%  Gaussian  noise  with  zero 
mean.  Figures  6  and  7  show  that  we  have  robust  convergence  as  well.  As  shown  in  Figure 
6(b),  we  are  not  able  to  drive  the  energy  to  zero.  This  is  expected  because  the  boundary 
measurements  are  highly  oscillatory,  and  in  general  we  cannot  find  a  smooth  velocity  c  which 
produces  exactly  the  same  travel-times  as  those  noisy  data. 

5  Three-Dimensional  Numerical  Examples 

In  the  following  examples,  we  use  65  x  65  x  65  grid  points  in  the  three-dimensional  space. 
Using  the  above  formulation,  we  need  measurements,  denoted  by  T*,  on  the  boundary  dflp. 
For  each  velocity  model  shown  below,  we  have  implemented  the  following  case,  49  sources 
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Figure  2:  (Example  2.  Ten  Sources)  (a):  the  relative  error  in  the  solution  and  (b):  the 
convergence  history  of  energy. 


on  the  levels  z  =  0.1  and  0  =  1.9,  and  we  have  98  sets  of  measurements  in  total. 

To  start  the  algorithm,  we  initialize  the  velocity  c°  by  solving  the  elliptic  equation  (18) 
with  v  —  1. 


5.1  Example  1.  Constant  Model. 

The  exact  velocity  model  is  given  by  c  =  1. 

We  use  the  BFGS  method  to  invert  for  the  velocity.  The  results  are  shown  in  Figure  8; 
we  observe  the  quadratic  convergence  once  again. 


5.2  Example  2.  Gaussian  Model. 

The  exact  velocity  model  is  given  by 


c(z,  y,z)=  3  -  -  exp 


x2  +  y2  +  (z  —  0.5)2 
052 


exp 


x2  +  y2  +  [z  —  1.25)' 
052 


(32) 


We  use  the  gradient  descent  method  to  invert  for  the  velocity.  The  results  are  shown  in 
Figure  9. 


6  Synthetic  Marmousi  model 

The  Marmousi  model  from  the  1996  INRIA  Workshop  on  Multi-arrival  Travel-times  is  a 
synthetic  model  which  will  challenge  the  adjoint  state  method  used  here. 
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(a)  (b) 

Figure  3:  (Example  2.  Ten  Sources)  Cross-sections  of  the  solutions,  (a):  z  —  1  and  (b): 
x  —  0. 

The  original  Marmousi  model  is  sampled  on  a  24m  by  24m  grid,  consisting  of  384  samples 
in  the  ^-direction  and  122  samples  in  the  ^-direction;  therefore  the  model  dimension  is 
9.192km  long  in  the  x-direction  and  2.904km  deep  in  the  ^-direction. 

In  the  computational  results  presented  here,  we  use  20  sources  and  their  (x,  ^-coordinates 
are  (200,  2800),  (1000:1000:9000,  2800),  (200,  100)  and  (1000:1000:9000,  100),  respectively, 
where  we  have  used  by  now  the  standard  Matlab  colon  notation. 

The  true  Marmousi  velocity  model  is  illustrated  in  Figure  10(a).  As  we  can  see,  this 
velocity  model  has  high  contrast  with  variations  of  different  scales.  On  the  one  hand,  since  the 
fast  sweeping  method  used  here  is  unconditionally  stable,  the  forward  eikonal  solver  will  not 
have  difficulty  in  computing  traveltime  to  first  order  accuracy.  On  the  other  hand,  viscosity- 
solution  based  first-arrival  traveltimes  will  not  be  able  to  give  us  too  much  information  about 
variations  of  small  scales  occurring  in  the  velocity  model;  to  retain  the  information  related  to 
small  scales,  we  have  to  use  multiple  arrivals,  which  in  turn  calls  for  multiple-arrival  based 
traveltime  tomography.  In  this  regard,  for  computing  multiple  arrivals  of  the  Marmousi 
model  in  the  Eulerian  framework,  see  [17,  18]  for  more. 

To  start  the  algorithm,  we  initialize  the  velocity  c°  by  solving  the  Laplace  equation 
—Ac0  =  0  with  c° |ao  =  cexactldn-  The  solution  is  plotted  in  Figure  10(b). 

We  use  the  BFGS  method  to  invert  for  the  velocity.  Figure  11  presents  the  inversion 
results  for  different  cases  in  terms  of  the  sampling  size  Ax  and  the  parameter  v. 

Comparing  Figure  11(a)  with  the  true  model  Figure  10(a),  we  have  succeeded  in  imaging 
the  macro  scale  variations  of  the  velocity  model  and  we  were  not  able  to  image  finer  scale 
variations  of  the  velocity  model  which  does  exist  in  the  true  velocity  model.  However, 
transmission  tomography  usually  has  very  limited  resolution,  and  we  believe  that  this  result 
is  near  optimal  using  the  current  approach. 

To  confirm  this,  we  refine  the  velocity  model  by  doubling  the  number  of  grid  points  in 
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(a)  (b) 

Figure  4:  (Example  3.  Ten  Sources)  BFGS.  (a):  the  final  approximated  c  and  (b):  the 
convergence  history  of  energy. 


each  direction  while  keeping  the  regularization  parameter  v  fixed;  the  corresponding  solution 
is  shown  in  Figure  11(d).  We  also  check  the  following  residual  in  the  solution  defined  by 


R 


N  fgn  ds 


(33) 


where  N  is  the  number  of  sources  defined  above.  This  quantity  essentially  is  the  average 
relative  error  in  the  first-arrival  time  per  source  per  receiver.  If  the  above  residual  is  not 
changing  too  much,  we  will  accept  the  inversion  result  since  there  is  not  much  model  misfit 
left  to  drive  improvement.  Indeed,  as  shown  in  Figure  12,  even  if  we  refine  the  velocity 
model,  the  residuals  are  almost  the  same  after  15  BFGS  steps.  In  fact,  the  solutions  from 
the  coarse  and  fine  resolution  are  similar,  as  shown  in  Figures  11(a)  and  11(d)  . 

Using  a  relative  small  v  =  102,  the  BFGS  iteration  has  difficulty  in  converging  to  a 
smooth  solution.  This  is  clearly  seen  in  Figure  11(b).  The  BFGS  iteration  stops  at  the 
fourth  iteration  with  E{mA)  ~  1700,  where  m  =  logc;  see  Figure  12.  This  difficulty  comes 
from  the  sharp  spikes  in  the  solution  near  the  source  locations,  where  the  traveltime  field  is 
not  differentiable  [19].  These  sudden  changes  will  degrade  accuracy  in  the  computed  gradient 
and  make  it  hard  for  BFGS  to  search  for  a  descent  direction. 

Increasing  the  magnitude  of  v  (from  104  to  106),  on  the  other  hand,  we  have  a  little  bit 
better  convergent  result.  As  seen  in  the  energy  plot,  the  energy  which  uses  the  larger  v  (the 
dashed  line)  reaches  a  lower  state  than  that  using  v  =  104  (the  solid  line).  Theoretically  we 
penalize  the  gradient  of  c  so  that  it  is  small  in  a  weighted  Sobolev  space  as  illustrated  in 
equation  (12). 

Concerning  the  speed,  the  computational  time  for  the  cases  with  v  =  104  using  Ax  =  24 
and  12  are  53  minutes  and  387  minutes,  respectively.  We  also  list  in  Table  1  the  number 
of  iterations  required  to  solve  both  the  eikonal  equation  and  the  adjoint  equation  for  each 
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Figure  5:  (Example  3.  Ten  Sources)  BFGS.  Cross-sections  of  the  solutions,  (a):  z  —  1  and 
(b):  x  —  0. 


Ax 

Eikonal  equation 

Adjoint  equation 

24 

12 

20  (6.68  x  10“1U) 
28  (9.60  x  10"11) 

17  (1.62  x  10~9) 
25  (3.38x“8) 

Table  1:  Iteration  count  for  the  fast  sweeping  methods.  The  numbers  in  the  brackets  are 
the  errors  in  the  corresponding  iteration,  ||Tn+1  —  Tn\  \  or  ||An+1  —  An||. 


given  velocity  held.  These  numbers  are  obtained  for  the  case  v  =  104  with  only  one  point 
source  located  at  (5000,  —2800)  in  the  first  BFGS  iteration.  The  first  row  shows  the  number 
of  iterations  with  Ax  =  24,  while  the  second  row  corresponds  to  the  case  with  Ax  =  12. 


7  Conclusion 

We  have  proposed  a  PDE-based  Eulerian  approach  to  traveltime  tomography  so  that  we 
can  avoid  the  cumbersome  ray-tracing  in  inversion.  We  started  from  the  eikonal  equation, 
defined  a  mismatching  functional  and  derived  the  gradient  of  the  nonlinear  functional  by  an 
adjoint  state  method.  The  resulting  forward  and  adjoint  problems  can  be  efficiently  solved 
by  using  the  fast  sweeping  method.  In  addition,  we  have  used  a  limited  memory  BFGS 
method  to  drive  the  functional  to  zero  with  quadratic  convergence.  Numerical  results  for 
2-D,  3-D,  and  Marmousi  synthetic  velocity  models  demonstrated  the  robustness  and  the 
accuracy  of  the  method. 

The  methodology  proposed  here  is  quite  general  and  can  be  extended  to  many  other 
situations  without  any  major  difficulty.  For  example,  instead  of  point  sources  we  can  easily 
modify  the  formulation  to  accommodate  plane  waves  as  the  source  condition,  which  can  be 
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Figure  6:  (Example  3  with  added  noise.  Ten  Sources)  BFGS.  (a):  the  final  approximated  c; 
(b):  the  convergence  history  of  energy. 


achieved  by  using  the  boundary  condition  resulting  from  the  plane  wave  condition  in  equation 
(2).  If  the  domain  to  be  imaged  is  irregular  or  non- rectangular,  then  we  can  use  the  fast 
sweeping  method  designed  in  [20]  to  solve  the  eikonal  equation  efficiently;  although  we  cannot 
directly  apply  the  standard  FFT  technique  when  regularizing  the  gradient  direction  (11), 
we  may  still  use  other  fast  solvers  like  multi-grid  methods  to  solve  the  Possion  equation. 
To  further  improve  the  resolution  of  inverted  velocity  models,  one  can  also  incorporate  the 
amplitude  information  into  the  formulation;  this  generalization  is  the  so-called  diffraction 
tomography,  which  consists  of  an  ongoing  project. 
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Figure  11:  (Marmousi  model)  Converged  solutions,  (a):  v  =  104  and  Ax  =  24;  (b):  v  =  102 
and  Ax  =  24;  (c):  v  =  106  and  Ax  =  24;  (d):  v  =  104  and  Ax  =  12. 
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Figure  12:  (Marmousi  model)  The  change  in  (I):  the  energy  and  (II):  the  residual.  Legend: 
(a):  v  —  104  and  Ax  =  24;  (b):  v  =  102  and  Ax  =  24;  (c):  u  —  106  and  Ax  =  24;  (d): 
v  =  104  and  Ax  =  12. 


